When aircrafts fly with angle of attack, since the pressure differential between the upper surface and the lower one lets air flows 6 rotate around the tip of wing 4, a continuous vortex generate at the wing-tip and go down stream, like shown in the FIG. 1. This flow pattern is called the wing-tip vortex flows 5, which are vortex-dominated flows. If aircraft flying speed is smaller than 0.3 Mach number, the flow is incompressible flow. Actually, the wing-tip vortex could damage aircraft stability and deduce the lift, because the wing-tip vortex can generate an unsteady downwash, which produces an induced drag to aircraft as to deduce the ratio of lift to drag, and more, which can spread down-stream as far as several kilometers and form a dangerous zone for the following aircrafts. Most aircrafts have the devices to reduce the wing-tip vortex.
There are several methods to study the wing-tip vortex, among which, the numerical simulation, combining fluid dynamics, applied mathematics, computer science, due to its low cost and visualization property, occupies an important role. It is an extended application of computational fluid dynamics (CFD) in exploration to the mechanism of fluid flows and in design of industrial products. A biggest challenge for CFD is to improve the accuracy of simulation, to reduce error, loyally descript the physical characteristics of the fluid flows.
One of the important factors to affect the accuracy of CFD is the numerical diffusion in the numerical solutions when numerical methods are used to solve the fluid flow governing equations, such as the Euler or the Navier-Stokes equations. For examples, the spatial discretization methods, such as central or upwind difference, the temporal discretization methods, such as explicit or implicit methods, the usage of turbulence models, such as the two-equation or LES model, and the orthogonality of computing grids, all could produce the numerical diffusion. Additionally, some artificial diffusion is often added in numerical solutions to improve the numerical convergence by deducing the computing accuracy. For examples, the capture of the discontinuity interfaces, such as shocks, is achieved by appropriately adding some artificial diffusion in avoiding the numerical oscillation around the large gradient in the numerical solutions with higher accuracy. The numerical diffusion can be considered as a kind of energy loss in flow field, which makes the computational results be not very accuracy. The Advanced numerical methods should minimize the numerical diffusion under premise of achieving the convergent numerical solutions.
The most significant effect of the numerical diffusion to numerical simulation results is on the capture to the discontinuity interfaces. As mentioned earlier, it physically makes sense to add some artificial diffusion to capture the strong continuity interfaces, such as shocks, since across which the velocity, density and pressure all are discontinuous and there exists the energy losing in the form of entropy increasing. However, over large artificial diffusions could make the strong continuity interfaces smeared, which reduces the prediction accuracy of the shock intensity and position. Another type discontinuity in flow field is the contact discontinuity. It is a weak discontinuity, across which density and velocity in tangent direction is discontinuous while pressure and velocity in normal direction is continuous. This kind contact discontinuity, as the most significant property of the vortex-dominated flows and widely exists in engineering applications, such as the wing-tip vortex flows, is more difficult to capture compared to the strong discontinuity in numerical solutions, since very little amount of numerical diffusion could make the contact discontinuity interfaces smeared to a point damaging the accuracy. This is why developing the numerical methods to simulate the vortex-dominated flows is a big challenge for CFD workers.
A method to improve the simulation accuracy of vortex-dominated flows, such as the wing-tip vortex flows, is that solving the flow governing equations by using very fine computing grid, which costs more computing time. Moreover, the computing error could accumulate with the increasing of the number of computing grid cells. Another method is that intensifying a variable describing vortex flow movement, vorticity, by artificially adding physics models when solving the flow governing equations. For examples, adding the point vortex model could artificially increase vorticity. Locally solving the vorticity equation could reduce the vorticity losing in its transportation. However, the usage of those methods is quite limited, since the point vortex model is only for some simple flow situations where the position of vortex is known. Solving the vorticity equation is more complicated than doing the Euler or Navier-Stokes equations themselves.
In the earlier of 20 century, a vorticity confinement (CV) method was created by Steinhoff[1] to improve the simulation accuracy of incompressible vortex-dominated flows. Its working principle is that a source term as a vorticity-related body force is added in the momentum equations of the flow governing equations as to the vorticity diffusion can be removed from the numerical diffusions, which keeps the contact discontinuity interfaces in vortex-dominated flows out of smearing, so that the vortex structure can be captured more accurately. The details of the vorticity confinement are as follows. Firstly, the governing equations for the incompressible, invicid flows include the continuity equation and momentum equations, which are
                                          ∇                          ·                              V                →                                              =          0                ,                            (        1        )                                                                                    ∂                                  V                  →                                                            ∂                t                                      +                                          (                                                      V                    →                                    ·                  ∇                                )                            ⁢                              V                →                                      +                                          1                ρ                            ⁢                              ∇                ρ                                              =                      B            →                          ,                            (        2        )            
where, {right arrow over (v)} is the velocity vector, {right arrow over (v)}=ui+vj+wk, including three-dimensional components u, V, w in the direction i, j, k in the Cartesian coordinate system; ρ, p and t respectively represents the density, pressure and time.
In above equations, the operator ∇ presents
      ∇          =                                    ∂                          ∂              x                                ⁢          i                +                                            ∂                                                                                  ∂              y                                ⁢          j                +                                            ∂                                                                                  ∂              z                                ⁢          k                      ;the symbol • does the inner product. In the vorticity confinement, {right arrow over (B)}, the source term in form of body force in the right hand side of the momentum equation (2), is defined as{right arrow over (B)}={right arrow over (n)}c×{right arrow over (ω)},  (3)
where, {right arrow over (ω)} presents the vorticity, {right arrow over (ω)}=ωxi+ωyj+ωzk in the Cartesian coordinate system; the operator x presents the cross product. Based on the definition of the vorticity,
            ω      →        =                  ∇                  ×                      V            →                              =                                                                            i                                            j                                            k                                                                                                          ∂                                                                                                                      ∂                    x                                                                                                                    ∂                                                                                                                      ∂                    y                                                                                                                    ∂                                                                                                                      ∂                    z                                                                                                      u                                            v                                            w                                                              =                                            (                                                                    ∂                    w                                                        ∂                    y                                                  -                                                      ∂                    v                                                        ∂                    z                                                              )                        ⁢            i                    +                                    (                                                                    ∂                    u                                                        ∂                    z                                                  -                                                      ∂                    w                                                        ∂                    x                                                              )                        ⁢            j                    +                                    (                                                                    ∂                    v                                                        ∂                    x                                                  -                                                      ∂                    u                                                        ∂                    y                                                              )                        ⁢            k                                ,where, {right arrow over (n)}c presents the gradient direction of the vorticity,
                                                        n              →                        c                    =                                    ∇              ϕ                                                                    ∇                ϕ                                                                  ,                            (        4        )            where, φ is the magnitude of the vorticity and |∇φ| is the magnitude of the gradient of the magnitude of vorticity. Then,
                              ϕ          =                                                                  ω                →                                                    =                                                            ω                  x                  2                                +                                  ω                  y                  2                                +                                  ω                  z                  2                                                                    ,                            (        5        )                                                                ∇            ϕ                                    =                                                                              (                                                            ∂                      ϕ                                                              ∂                      x                                                        )                                2                            +                                                (                                                            ∂                      ϕ                                                              ∂                      y                                                        )                                2                            +                                                (                                                            ∂                      ϕ                                                              ∂                      z                                                        )                                2                                              .                                    (        6        )            
The source term {right arrow over (B)} in the momentum equation (2) needs to be multiplied by a factor ε,{right arrow over (B)}=ε({right arrow over (n)}c×{right arrow over (ω)}),  (7)which physically means a force at the center point of the vortex. In Steinhoff's original vorticity confinement, the factor ε is a constant value (0.01˜0.05), which can adjust the level of the vorticity confinement. This method has a simple formulation and needs not fined computing grids, so that it had been widely used in the simulation of vortex-dominated flows. The modification to this method in the following years is focused on the expression of the factor ε[2][3]. Although the vorticity confinement, based on the principle of adding a certain amount of vorticity to counteract the numerical diffusion on the contact discontinuity, had significantly improved the accuracy of simulation in vortex-dominated flows, it has the following defects,
1. The added vorticity is totally controlled by the factor ε;
2. The order of the spatial discretization of the added vorticity is limited;
3. The effect of the added source term on the convergence of the momentum equation is uncertain.
To more accurately simulate the wing-tip vortex flows, it is necessary to modify the mechanism to capture the contact discontinuity in incompressible vortex-dominated flows field. A new way is to refine the precision of the added vorticity in the source term in the momentum equation to counteract the numerical diffusions. This invention presents a new numerical method, which is called Vorticity-Refinement, where undertakes the high-order spatial discretization to the added vorticity and improve the convergence effect of the source term on the momentum equation.